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\0 : ABSTRACT 
O ■ 

Qh| The formation of vortices in protoplanetary disks is explored via pseudo- 

O ' spectral numerical simulations of an anelastic-gas model. This model is a cou- 

^ ' pled set of equations for vorticity and temperature in two dimensions which 

■ includes baroclinic vorticity production and radiative cooling. Vortex formation 

■ is unambiguously shown to be caused by baroclinicity because (1) these simula- 
^ ' tions have zero initial perturbation vorticity and a nonzero initial temperature 

■ distribution; and (2) turning off the baroclinic term halts vortex formation, as 
shown by an immediate drop in kinetic energy and vorticity. Vortex strength 
increases with: larger background temperature gradients; warmer background 
temperatures; larger initial temperature perturbations; higher Reynolds number; 
and higher resolution. In the simulations presented here vortices form when the 
background temperatures are ~ 200K and vary radially as r~°'^^, the initial vor- 
ticity perturbations are zero, the initial temperature perturbations are 5% of the 
background, and the Reynolds number is 10^. A sensitivity study consisting of 
74 simulations showed that as resolution and Reynolds number increase, vortices 
can form with smaller initial temperature perturbations, lower background tem- 
peratures, and smaller background temperature gradients. For the parameter 
ranges of these simulations, the disk is shown to be convectively stable by the 
Solberg-H0iland criteria. 
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Subject headings: accretion, accretion disks, circumstellar matter, hydrodynam- 
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Introduction 



Planetary formation theories come in two flavors: core accretion, where particles build 
from millime ter-sized dust g rains to kilometer-sized planetesimals and eventually to plan- 
etary cores (jWetherill 1990 ): and gravitational instability, where part of a mass ive disk 
quickly collapses to form giant planets by gravitational self-attraction (jBossI 119971 ) . Each 
theory has its own difficulties. Core accretion is an extremely slow process — it is thought 
that core accretion wou ld take 8 to 48 million years to build a Jupiter-sized planetary 
core (jPoUack et al.l 119961). This is i nconsistent with the estimated disk life-time of a few 
million years (IBricefio. et. al.ll200ll : iBally et al.lll998l ). Also, the mechanism for the cohe- 
sion of meter-sized particles is poorly understood; they are too small to be influenced by 
gravity, but too large for agglomeration by collision al coagulation in a turbulent gas disk 
( 1 Weidenschilling] 1 1984 IWeidenschilling fc Cuzzilll993l ). Gravitational instabilities could form 
planets fast enough, but the densities required for instabilities are hard to attain: high gas 
densities throughout a thick disk would produce global gravitational instabilities that would 
transport a substantial fraction of the mass into the central star; a thin, high density layer 
of solids in the rn idplane of the disk may be precluded by shear turbulence in the vertical 
dCuzzi et al.lll993h . 



Protoplanetary disks initially consist of about 99% gas and 1% dust ( IHayashil Il98ll ). 
Although most discussions of planetary formation only consider solid materials, the gaseous 
dynamics can have a large influence on the distribution of the solids and thus the subse- 
quent pla net building. Vorti c es in the gas disk are extremely efficient at capturin g meter-size 
particles ( Tanga et al. 1996 : Johansen et al. 2004 : Klahr fc Bodenheimer|]2006l ). so particle 
concentrations inside vortices would be several orders of magnitude higher than in the sur- 
rounding gas (IBarge Sz Sommerialll995l ). Thus vortices would be beneficial for either theory: 
in core accretio n, the rate of accret ion depends strongly of the local concentration of parti- 
cles in the disk (jPoUack et al.lll996l ): similarly, a high local density of solids within a vortex 
could initiate a gravitational instability without requiring a high density throughout the disk. 
The disks considered here are cool enough that they are only weakly ionized, and therefore 
turbulence due to magnetohydrodynamic effects is not considered. 

In order to know whether vortices play a role in planetary formation, we must under- 
stand what conditions are required for vortex formation a nd longevity. Several such models 
have been introduced in the literature in the last decade. iBracco et al.l (119991 ) used a two- 
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dimensional incompressible shallow water model with a background Keplerian rotation to 
study the inverse energy cascade. They found that if the initial energy of the disturbance 
is less than about 10~^ of the energy in the Keplerian flow, the vorticity fluctuations shear 
away and the disk returns to its unperturbed velocity profile. For lar ger initial energies, an - 
ticyclonic vortices form and merge, while cyclonic vortices shear out. iGodon fc Livid (119991 ) 
studied vortex longevity with a compressible Navier-Stokes pseudospectral model. Vortex 
life-time, as measured by the maximum vorticity, was highly dependent on the Reynolds 
number; at Re = 10^-10^ the maximum vorticity had an e-folding time of 60 orbital periods. 
They also found that Keplerian shear inhibits the formation of vortices larger than the scale 
length Ls = \J v / [dj.^) , where v is the rotational velocity of the vortex and Vt is the angular 
velocity of the disk. 

The only source of vorticity in two-dimensional (2D) models is baroclinicity, which 
appears as Vp x Vp in the vorticity equation, where p is the pressure and p is the density. 
In the atmosphere, baroclinicity usually operates in the vertical and produces effects such 
as land-sea breezes when density and pressure (or temperature) gradients do not align. In 
the geometry of a protoplanetary disk, density stratification is in the radial, so azimuthal 
variations in pressure will produce vorticity. Most numerical models of protoplanetary disks 
do not include the effects of baroclinicity and can therefore only simulate decaying turbulence 
fro m some initia l vorti city distribution. For ex a mple, the incompressible fluid models used 



bvlBracco et all f|l999l ) and 
Godon fc LiJiTjl999h and 

p = Kp^ so that Vp X Vp = 



Umurhan fc Reeev 



2004r) p reclude any thermodynamic effects; 



Johnson fc Gammid (120051 ) both assume a polytropic relation 



and consequently baroclinic effects are not present. 



Vortex formation due to baroclinic effects was studied by lKlahr fc Bodenheimerl (120031 ). 
who found that the generation of turbulence and vortices depends strongly on the radial 
temperature gradient. They observed vortices in both three-dimensional (3D) simulations 
with a full radiative scheme (a flux-limited radiative transport with a simple dust opacity), 
and 2D simulations where the radiative transport scheme was replaced by a specifled back- 
ground temperature distribution. When the radial temperature was constant, no turbulence 
w as produced; when T ^ r~^, t he flow became turbulent and vortices formed. The results 
of iKlahr fc Bodenheimerl (12003! ) are the motivation for the current study; we aim to eluci- 
date how baroclinic vorticity production depends on background temperature, background 
temperature gradient, and the thermal diffusion. Both our work and the simulations in 
Klahr fc Bodenheimerl (120031 ) have disks which are convectively stable in the radial, based on 
the Solberg-H0iland criteria. 



The study of 3D vortices by iBarranco fc MarcusI (120051 ) models an anelastic fluid using 
the Euler equation and the ideal gas law; thus baroclinic effects are included. They found 



-4- 



that vortices that extend several scale heights vertically are unstable, as are short vortices 
centered at the midplane. However, breaking internal gravity waves would spontaneously 
create 3D vortices centered 1-3 scale heights above or below the mid-plane. This has obvious 
implications for the stability of the vortices that we observe in our 2D model, which are 
similar to their tall unstable vortices. This issue is discussed further in the conclusion. 



Notably, even though iBarranco fc MarcusI (120051 ) include baroclinic effects, they do not have 
a radial temper ature gradient in their mode l, so their vortices are not produced by the 



Klahr fc Bodenheimerl (120031 ) or ours, which both require a nonzero radial 



same process as 
temperature gradient. 

This work is presented in two parts, where Part I studies vortex formation with small 
domain simulations for five orbital periods, and Part II studies vortex growth and longevity 
using a larger domain for hundreds of orbital periods. This paper begins by describing 
the 2D anelastic gas equation set in ^ and the pseudo-spectral numerical model in ^ The 
results in §1] include a description of the early evolution of these simulations, and a sensitivity 
analysis of how vortex formation depends on numerous parameters. Conclusions follow in 
§5] The appendix includes details of the numerical model, including enforcement of stress- 
free boundary conditions with the influence matrix technique. For conciseness there is little 
repetition between Parts I and II, so the reader is advised to read both together. 



2. Description of the Equation Set 

Any numerical model of an astronomical object must make compromises between the 
competing demands of including the essential physical processes, achieving adequate nu- 
merical resolution, and exploring the relevant range of parameters. Given the con troversy 



generated by the large scale simulations published by iKlahr fc Bodenheimerl (120031 ). we be- 
lieve it is prudent to use a simplified model that can help elucidate the necessary requirements 
for baroclinic vorticity production in the 2D disk geometry. More realistic models can be 
used to verify the relevance of our model results to actual protoplanetary disks. 

In this study we model a 2D anelastic gas, rather than compressible gas, in order to 
make a fast numerical model that can explore a large swath of parameter space by performing 
numerous simulations. Anelastic models filter out the acoustic waves present in compressible 
gases; the lack of these fast-propagating acoustic waves relaxes the time step restrictions of 
the numerical modes, thereby reducing the computational cost. By using a fast numerical 
model, we are also able to run simulations at higher resolutions than most published models 
of astrophysical disks and thereby minimize the effects of numerical viscosity. Compressible 
disk simulations typically have such large numerical viscosities that they can suppress the 
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baroclinic production of vorticity we want to explore (IKlahr fc Bodenheimerll2003l ). A fully 
3D anelastic model would also be very expensive to run at the same resolution that we can 
easily explore with our 2D model. Finally, our model retains the essential physics of interest, 
namely the baroclinic instability and radiative cooling. 

To arrive at the anelastic limit, assume that thermodynamic perturbations are small 
compared to the background state and that radial displacements are significantly less that 
the le ngth scale for radial gradients in the unperturbed surface density of the disk (IBannon 
19961 ). Then the conservation of mass equation. 



Dp 



u 



0, 



simplifies to the anelastic equation. 



V ■ (pu) = 0, 



(2) 



where the density p is constant in time and varies only in the direction of stratification. Here 
u is the velocity vector and D/Dt is the material derivative. 



The model equation set is inspired by iBannonI fjl996l ) and lScinocca fc Shepherd! (119921 ). 
which are anelastic models of the atmosphere derived from conservation of momentum, 
conservation of mass, the second law of thermodynamics, and the ideal gas law. These 
equations were modified to the geometry of a disk, where the background density and tem- 
perature stratification are in the radial and shear due to near-Keplerian rotation is in the 
azimuthal. Assuming that the disk is thin and hydrostatic in the vertical leads to a 2D model. 
A vorticity-stream function formulation is used so that only two prognostic variables, the 
vorticity and the potential temperature, must be computed by the numerical model. 



The model equations are 



c 



09 
dt 



Id r r 
r dr \So dr 
Cp diTo 86' 
r dr d(f) 
9' 

T 



1 



r2Sr 



(3) 
(4) 
(5) 



The vertical component of vorticity, the stream function, and the potential temperature are 
the sum of a background and perturbation variables, 

C = Co(r) + C'(r,0,t) 

^ = *o(^) + *'('^,0,^) 

9 = 9o{r) + 9\r,(P,t), 
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where the primed perturbation quantities are small relative to the background. Radial and 
azimuthal velocities u = [u, v) are related to the streamfunction by Squ = —V x ^z. Other 
variables include the surface density So(r), Exner pressure 7ro(r), radiative cooling time r, 
specific heat at constant pressure Cp, time t, polar coordinates (r, 0), vertical unit vector 
z, and the Jacobian d{a,b) = {drad^h — d^adrb)/r. The potential temperature and Exner 
pressure can be written in terms of the pressure p and temperature T as 

\ P J d \Po{rin)J 

where po(^m) is a reference pressure, rj„ is the inner radius of the annulus, and R is the gas 
constant. The background potential temperature profile, 9o{r), is assumed to be maintained 
by a balance between radiative heating by the central star and radiative cooling to space. 

The background variables obey a static background balance between the centrifugal 
force, gravity, and the pressure gradient, 

where Qq is the background angular velocity, $ is the gravitational potential, and po is the 
pressure integrated over the thickness of the disk. This leading order balance has already 
been subtracted out of the vorticity equation (j4]). 

Equation (jl]) is a conservation equation for potential vorticity, C/^O) where the baroclinic 
term on the right-hand side is a source of vorticity and the cause of the baroclinic instability 
that is central to our study. This baroclinic term will look unfamiliar to most readers due 
to the variables and geometry used. Recall that the common form comes from taking the 
cross-product of the pressure gradient term in the momentum equation, 

V X [--Vp] = X Vp. (7) 

V P / P 

The equivalent operation in our case is 

I d f ^,d'^o\ Cpd9' drvQ 



r d(j) dr 



In the numerical model, effective eddy dissipation terms are added to the vorticity and 
temperature equations for numerical stability, and many of the background variables drop 
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out: 



^, _ 1 (9 / r d^'\ ^ 1 92^' 
r dr \TiQ dr J r^Eo 



|.,i_a,*,„ . 4 + (11) 

Here z/g is the viscosity and Kg is the thermal dissipation. 

The equation set is then nondimensionahzed using the length scale = rmid'-, the time 
scale = 27r/f2o (^^mid); which is one orbital period; and typical temperature and density 
scales. For the standard domain, Lgc = 7.5AU and tsc = 20 (earth) years, background 
temperatures are 125-600K, and background surface densities are 350-1000g/cm~^. These 
scales let us define the Reynolds and Peclet numbers based on the eddy dissipation terms, 

7-2 7-2 

Re = Pe - 



which measure the ratio of momentum advection to kinetic energy dissipation and thermal 
advection to thermal dissipation, respectively. For a particular grid spacing. Re and Pe 
must be small enough to avoid building up energy at the smallest scales. Simulations with 
a larger number of gridpoints allow larger values of Re and Pe, which reduces dissipation so 
that coherent structures are more long-lived. For brevity, the nondimensionahzed equations 
are not presented here. 

Factors of the background surface density Eg appear in the model equations because \1/ 
is a momentum stream function, defined as Equ = —V x where u = {u,v), the radial 
and azimuthal velocities, and z is a vertical unit vector. A momentum stream function 
was chosen because the anelastic equation, V ■ (Equ) = dictates that the product EqU 
is divergence free. (In incompressible fiows, V ■ u = and the stream function is simply 
u = —V X \l/z.) Equation ([9]) is nearly a polar Laplacian (' = V^\l/'. The density factors 
present an extra challenge for the inversion step to find (appendix IA.3p . The Jacobians on 
the left-hand side of (fTO!) and (ITT!) represent a source in the Eulerian frame due to advection, 
and combine with the time derivatives to form material derivatives in polar coordinates. 

Equation fill I) is an advection-diffusion equation for potential temperature perturbation, 
6'. The two terms on the right-hand side are both decay terms: the Laplacian operator 
imparts an effective eddy diffusion which dissipates more quickly at the smallest scales; the 
radiative damping term, —6'/t, is a simplified model of blackbody radiation and dissipates 
perturbations equally at all scales so as to relax the temperature field back to the background 
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value. The background potential temperature profile, 6o{r), is assumed to be maintained by 
a balance between solar heating and radiative cooling to space. 

The temperature equation flTTl) is coupled to the vorticity only by its advection (Ja- 
cobian) term, while the vorticity equation (fTOl) is coupled to temperature by its baroclinic 
term. This coupling allows the baroclinic feedback between vorticity and temperature per- 
turbations to occur, which can cause nonlinear growth leading to turbulence. Understanding 
what parameters affect the strength of this baroclinic feedback is one of the main goals of 
this work. 



3. The Numerical Model 



A two-dimensional spectral model was created to investigate the nonlinear behavior of 
the equation set. Spectral methods have the advantages of exponentially convergent spatial 
representation and efficient inversion of derivative operators. Vorticity and temperature fields 
on the annulus are represented using Fourier basis functions in the azimuthal and Chebyshev 
basis functions in the radial. The time-stepping algorithm is a third-order semi- implicit 
Runge-Kutta routine which is optimized for efficient memory use ( ISpalart et al.lll99ll ). Other 
authors that have cre ated spectral models on an annulus o r dis k using Fourier- Chebyshe v 
basis functions include iGodonI (119971 ) . iBergeron et al.l (120001 ) . and lTorres fc CoutsiasI (119991 ). 



The Laplacian operator, which is inverted at each stage of the semi-implicit Runge-Kutta 
routine, is an upper triangular matrix when applied to the Chebyshev coefficients for each 
Fourier mode. This upper triangular matrix was reduced algebraically to a banded diagonal 
of width six, greatly improving the efficiency of the inversion step, which is potentially a very 
time consuming operation (see appendix I A. 2|) . Specifically, inverting the Laplacian operator 
with N Chebyshev modes and M Fourier modes, the operation count is 0{MN'^) for a 
back-solve routine on a triangular matrix, and is reduced to 0{MN) for a banded diagonal. 

A standard polar Laplacian operator is inverted for the vorticity and temperature equa- 
tions, (fTOj) and (ITT!) . For the inversion of the near-Laplacian with surface density factors in 
([9]), it is not possible to solve for using an arbitrary function So(r). However, a simple 



power function Sq 



cr 



only adds a constant in front of one of the radial derivative terms. 



as described in the appendix IA.3[ This method allowed us to incorporate radial variations 
in the surface density. The restricti on to power laws is not overly prohibitive, as they are 
standard test cases in other studies ( Klahr fc Bodenheimer 2003 : Godon fc Liviol 1999 ). 



All derivatives are computed in spectral space using the Fourier- Chebyshev coefficients. 
The products that appear in the advective and baroclinic terms are computed at each grid 
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point in physical space. In this pseudospectral method, the vorticity and temperature fields 
must be transformed between p hysical and spectral sp ace at every step. We use the FFTW 



fast Fourier transform package flFrigo fc Johnsonlll998l ). which has been benchmarked as one 



of the fastest FFT algorithms on numerous platforms. 

The first task in running fiuid simulations is to find the appropriate dissipation rate for 
the grid spacing. The goal is to have as little dissipation as possible (i.e., large Reynolds and 
Peclet numbers), while dissipating enough energy at the smallest scales to prevent grid-scale 
oscillations and numerical instabilities. Calculating a one-dimensional energy spectrum in 
Fourier- Chebyshev space is not straightforward, so separate plots of energy versus Fourier 
modes and Chebyshev modes were checked to ensure that the energy decays properly at 
small scales. 

Additional viscosity was added at large radial wavenumbers to dissipate energy near the 
gridscale, so that the radial Laplacian terms in (fTOj) and (fTT!) have the coefficient 



1 

''coef 



2 



Uf + l + {uf- 1) tanh {^-^^ 



(12) 



in spectral space. For wavenumbers somewhat smaller than the critical wavenumber fcc, 
i^coef ~ 1 and the diffusion operators are unaffected; for wavenumbers k » kc the coefficients 
1/Re and 1 / Pe are multiplied by the factor i^coe/ ~ J^/- Typical simulations used a factor Uf oi 
500, a critical wavenumber kc of four-fifths the maximum wavenumber, and a transition width 
kw of 10. These parameters were chosen so that only the largest Chebyshev (radial) basis 
functions experience additional damping. It was found that hyperviscosity was not needed 
for the Fourier (azimuthal) basis functions. With the addition of this hyperviscosity, we were 
able to increase the effective Reynolds number from 5 x 10^ to 4 x 10^, and develop much 
finer-scale structures at the same resolution. Including hyperviscosity does not affect the 
conclusions of this study; all the simulations in Tabled] were also run without hyperviscosity 
and produced the same qualitative results. 

The azimuthal domain can be chosen from a full annulus or any fractional annulus 
27r/2", such as a quarter or 64th. Fourier basis functions are periodic in the azimuthal in all 
of these cases, but the indexing of the basis functions and their derivatives must be handled 
with care. For example, if a function is represented on the full annulus as 

N 



Y^aue'^^ (13) 



fe=0 



then on the half annulus all of the odd modes must be zero, so only the even modes 
Oo, a2, 04, ag . . . survive. For efficient memory storage, only the nonzero coefficients are kept 



-lo- 



in the array, so the multiphcation factor of a derivative d/d(f) ^ ik is not the same as the array 
index. The simulations in this paper are on a 64th annulus, with a resolution of 128 x 128 or 
256 X 256. Paper II presents simulations on the quarter annulus with resolutions of 256 x 256 
and 512 x 512. 



3.1. Boundary conditions 

Fourier basis functions are the natural choice for the periodic boundary conditions in 
the azimuthal direction of the annulus. The radial perturbation temperature boundary con- 
dition is zero flux (Neumann), drO' = 0, at both the inner and outer boundaries. Physically, 
the zero flux conditions means that the disk is not influenced by heating from outside of 
the computational domain. They are enforced during the inversion of the Laplacian oper- 
ator using the tau method, where two rows of the Laplacian matrix are substituted with 
the boundary constraints. This effectively means that the projection of the temperature 



perturbation variable onto the highest t wo Chebyshev 



constraints on the boundary conditions (jFornberglll996l ) 



3olynomial modes are exchanged for 



The dynamic boundary conditions are stress-free and preclude flow normal to the inner 
and outer boundary, 

'^-- = 0, „' = 0. (14) 

or r 

This choice is sensible for astrophysical applications, as there are no solid surfaces to impart 
stress at the boundaries. Inflow/outflow conditions that drive the perturbation flow through 
boundary forcing are not considered in this study. The stress-free boundary conditions are 
rewritten in terms of the stream function and vorticity, and are enforced using the influence 
matrix technique (appendix IA.5p . 



3.2. Initial Conditions 

Initial conditions used for the perturbation variables (' and 6' include: (1) Fourier- 
Chebyshev modes; (2) Fourier- Bessel modes; (3) Fourier-sine modes; (4) a collection of 
vortices; and (5) a specifled energy spectrum which is random in phase. Individual modes 
were used extensively in testing the code, as described in the model validation section 13.31 
Shearing effects and orbital time periods are clearly visible with the Fourier-sine modes. 
Individual vortices were used to compare the behavior of cyclonic versus anticyclonic vortices, 
and to observe vortex merger. 
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(Mc Williams 


1990; 


Bracco et al. 


2000) 



model was used to create these initial vorticity or temperature fields. The initial energy 
spectrum is E{k) ~ /c"/(fc^" + /cg")' where k is the Cartesian wave number modulus and ko 
and a are parameters that change the shape of the spectrum. The simulations discussed in 
this paper use ko = 10 and a = 10, as shown in the initial temperature perturbation in Fig. 
m Sensitivity studies with showed that the qualitative characteristics of the simulations did 
not depend on the particular value of ko. 

High resolution initial conditions from the Cartesian model were interpolated to the 
Fourier- Chebyshev gridpoints of the annulus using bilinear interpolation. The inner and 
outer radial edges were interpolated to zero to conform with boundary conditions, and the 
azimuthal edge was linearly interpolated to a periodic azimuthal boundary. 



3.3. Model validation 

The numerical model was validated by direct comparison with exact linear solutions 
using Fourier-Bessel functions. In the absence of the nonlinear and baroclinic terms, ffTOl) 
and flTTl) become heat equations in vorticity and temperature. For a single Fourier-Bessel 
basis function of the same mode number, 

V'Mr)e''"^ = -4(r)e*'=^ (15) 

where Jk{r) is the A;*^ order Bessel function of the first kind in the radial direction and 
^ikip ^YiQ f^th pQij^rigr mode in the azimuthal. Thus the exact solutions of (fTOj) and (fTTl) 
decay exponentially with e-folding times of 1/Re and 1/r + 1/Pe, respectively. Including 
the baroclinic term adds a non-homogeneous term to the exact solution of the temperature. 

As there are no exact solutions to the full nonlinear model, inclusion of the advection 
terms were tested in other ways. The Jacobian subroutine was tested individually and com- 
pared to analytic expressions to verify proper operation. When the advection is computed 
using only the background temperature, vorticity, and stream function, it becomes a linear 
term and an exact solution is available for comparison. Finally, qualitative results of vor- 
tex co-rotation, merger and translation showed that the numerical model was working as 
expected. 
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4. Results 

Hundreds of simulations were performed to study the effects of varying the model pa- 
rameters. The simulations discussed in this paper, shown in Table [1], vary parameters that 
affect vortex formation. The domain of these simulations have a radial and azimuthal extent 
of r G [9.5, 10]AU and (j) G [0,7r/32], respectively, and last for five orbital periods (here 
an orbital period is a full orbit of 2ti at mid-annulus, rmid = 9.75AU). Paper 11 discusses 
simulations with larger domains and longer duration, up to 600 orbital periods. 

The background surface density and temperature are power functions, 

So(r) = af— y, To{r)=c(—)\ (16) 

where is the inner radius of the annulus. All simulations presented here use a = 500g 
cm~^ and b = —3/2, while c and d are varied and shown in Table [H The resulting back- 
ground temperature distributions and corresponding potential temperature distributions are 
shown in Fig. [H All the remaining results in this paper are presented in terms of thermal 
temperature T, rather than potential temperature 6. This is a choice of convention; the 
model equations ([5]) are in terms of potential temperature, but thermal temperature is more 
intuitive and easier to compare to observations. 

The Schwarzschild criterion states that a disk is convectively sta ble in the absence o f 



rotation and shear if entropy, S, increases radially, i.e. dS/dr > (ISchwarzschildl Il958l ). 
Entropy is related to potential temperature as 

S = c,\n(^-^^ (17) 

where 6^ is a reference temperature, so the Schwarzschild criterion in terms of potential 
temperature is simply dd/dr > 0. The Schwarzschild criterion is satisfied for simulations 
where d > —1 (Fig. [T]). 

In the presence of differential rotation, a sufficient condition for convective stability is 
the Solberg-H0iland criterion, 

1 1 

-4 Vp-V5>0, (18) 

r'^ or CpP 



where j = r'^Q is angular momentum per unit mass (jTassou]|l2000l : iRiidiger et al.ll2002l ). In 
our variables, this criterion is 



dr 6'oSo dr dr 
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Differential rotation has a stabilizing effect on disks. In other words, disks which do not 
satisfy the Schwarzschild criterion may still satisfy the Solberg-H0iland criterion if the shear 
[dj'^/dr) is large enough. This is true for the simulations presented here; when d = —1 and 
d = —2, the left-hand side of f[T^ ranges between 0.034 and 0.044 years"^. In summary, all 
of the simulations discussed in this paper are convectively stable. 

The initial condition for the perturbation temperature was a specified distribution in 
spectral space (Sect 13.21) and a magnitude in physical space which is 5% of the maximum 
background temperature. The initial vorticity is zero; the only way for vorticity to be pro- 
duced in these simulations is through the baroclinic term. In a typical simulation, the initial 
temperature perturbations create vorticity perturbations in this way during the first several 
orbital periods (Fig. [2]). The maximum perturbation vorticity \('\, the perturbation kinetic 



energy + f'^|/2, and the maximum perturbation temperature \T'\ for simulation A are 
shown in Fig. [5l The vorticity and kinetic energy grow at early times due to the baro- 
clinic term, while the temperature perturbations die off due to thermal dissipation. After 
about one orbital period, the baroclinic vorticity production stops because the temperature 
perturbations are so low that dO' /d(j) in the baroclinic term is also small. During the first 
orbital period, differential rotation causes the vorticity perturbations to be stretched out 
into narrow bands of positive and negative vorticity. At the same time, the local radial 
gradient in the vorticity perturbation rapidly grows to a magnitude that is comparable to 
the vorticity gradient of the background Keplerian flow. This kind of shear-induced ampli- 
fication of vo rticity perturbatioiis in astrophysical disks has been described in many recent 
publi cations f Chagelishvili et al. 2003[ ^ Teyzadze et al. 2003; Yecko 2004; Umurhan fc Regev 



2004h . In particular, lAfshordi et aP J2005h show that this transient amplification occurs m 
three-dimensional disks so long as the radial wave number is large compared to both the 
azimuthal and vertical wave numbers. 

Once the radial profile of the total vorticity in the simula tion develops inflectio n points, 
it will satisfy the requirements for a secondary instability (IDrazin fc ReidI Il98ll ). In the 
simulation, this secondary instability manifests itself after the flrst orbital period by causing 
regions of negative vorticity adjacient to strong vorticity gradients patches to roll up into 
anticyclonic vortices (vortices rotating in the opposite direction as the background flow), 
while the positive vorticity regions tend to shear out and eve n tually diffuse away (Fig. [3]). 
Similar secondary instabilities have been reported by iLi et al.l (120051 ) when simulating disks 
perturbed by an embedded planet. The persistence of anticyclone s and the absence of 



cyclones is typical of simulat i ons w ith background differential rotation (IGodon &: Lividll999 



Marcus! Il990l : iMarcus et al.l l2000l ). The vortices are closely coupled to the temperature 
perturbations. Figure H] shows that each vortex advects cold fluid from the outer annulus 
inwards and advects warm fluid from the inner annulus outwards. This net radial transport 
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of heat down the background temperature gradient plays an important roile in the barochnic 
feedback and long-term evolution of anticyclones, as discussed in Part II. 

The baroclinic term is the only source of perturbation vorticity in the vorticity equation 
( ITO|) . Since these simulations begin with zero vorticity perturbation, the baroclinic term 
must be responsible for the increase in vorticity at very early times. To explore this further, 
simulation B repeated simulation A, but turned off the baroclinic term at t = 0.1,0.2,0.4, 
and 0.8 orbital periods. As soon as the baroclinic term is turned off, the vorticity and kinetic 
energy immediately drops (Fig. [5]). If the baroclinic term is turned off sufficiently early, no 
vortices form; the longer it stays on during the first orbital period, the stronger the vortices 
become. These results show that the baroclinic term is responsible for the production of 
vorticity at early times. 

The strength of vortex formation was found to depend on the following parameters: the 
background temperature gradient d, the magnitude of the background temperature c, the 
size of the initial temperature perturbation, the Reynolds number Re, and the resolution. In 
order to explore this large parameter space, a sensitivity study consisting of 74 simulations 
was conducted. These are simulations P, TG, and T in Tabled and the results are shown in 
Figs. [6l U\ and [H respectively. The contour plots of minimum vorticity show the strength of 
the vortices after five orbital periods. When all other parameters remain constant, stronger 
vortices result from higher Reynolds numbers and higher resolution. These factors are not 
independent — higher resolution is required for higher Reynolds number. At resolutions of 
128 X 128 and 256 x 256 the highest attainable Reynolds numbers were 10® and 10^, respec- 
tively. These are still far from realistic Reynolds numbers for a protoplanetary disk, which 
would be 10^^ or highei0. 

The contours in Fig. slant from the top left to the bottom right of the contour 
plot. This means that stronger vortices may be produced with either a stronger initial 
temperature perturbation or a higher Reynolds number. Stated another way, as Reynolds 
number increases (i.e. becomes more realistic) a smaller initial perturbation is required to 
kick off vortex formation. Thus even though we can produce vortices with a 5% initial 
perturbation at 256 x 256, the trend suggests that a much smaller perturbation would be 
required at higher resolutions and higher Reynolds numbers. The contour plots in Figs. [7] 
and [8] show similar trends for the background temperature gradient d and the magnitude 
of the background temperature c. As shown by the snapshots of vorticity perturbation. 



^ Goldreich fc WardI (Il973h estimate the kinematic viscosity of a protoplanetary disk at lAU to be t/ = 

2 X 10^ cm^s"-'^. Using a length scale of lAU and a time scale of one year, the Reynolds number Re = 
lyitiy) = 4 X 1012. 
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when the Reynolds number is sufficiently high (Re = 10^), vortices form in these simulations 
when To = 300K{r / 9.5 AU)'^-^^ (Fig. d point D) and when Tq = 200A'(r/9.5Af/)-°-5 (Fig. 
El point D). In summary, vortex strength increases with: larger background temperature 
gradients (more negative d); warmer background temperatures c; larger initial temperature 
perturbations; higher Reynolds number; and higher resolution. 

We can summarize the vorticity evolution displayed by our simulations with the follow- 
ing steps: (1) the initial potential temperature perturbations lead to baroclinic production 
of vorticity perturbations; (2) the differential rotation of the background flow causes the 
vorticity perturbations to be stretched out in the azimuthal direction and to develop large 
radial gradients; (3) the radial gradients in vorticity perturbations become large enough to 
trigger a secondary instability that is caused by local inflection points in the radial pro- 
file of the total vorticity. The secondary instability causes regions of negative perturbation 
vorticity to roll up into nonlinear anticyclones. A major result of this paper is our finding 
that insufficient numerical resolution or too large a numerical viscosity will cause the vortic- 
ity perturbations to be damped out before step (3) can occur. Conversely, our simulations 
show that increasing both the numerical resolution and the effective Reynolds number of our 
simulation allows weaker initial temperature perturbations to evolve into nonlinear vortices. 

In Part II the radiative cooling time r and the Peclet number Pe are shown to play an 
important role in the baroclinic feedback. These parameters were varied in simulations Tau 
and Pe presented here to see their effects at early times. The variables r and Pe control 
the rate of dissipation in the potential temperature equation flTT]) : small values of r and Pe 
(fast diffusion) return the perturbation potential temperature 6' to zero faster than large 
values. This is exactly what is observed when either r or Pe are varied (Figs. [9] and [TOj) . 
As the temperature perturbations drop off, temperature gradients become smaller, and the 
baroclinic term produces less vorticity. Thus at early times, faster dissipation (smaller r or 
Pe) results in lower perturbation temperature, vorticity, and kinetic energy. At early times, 
vortex formation is strongest in the limit of no thermal dissipation (r, Pe oo). Part II will 
show that faster dissipation has the opposite effect during the vortex growth phase; smaller 
values of r or Pe result in stronger vortices due to the baroclinic feedback. 

5. Conclusions 

In this paper we have: introduced a two-dimensional anelastic equation set which is a 
simplified model of a protoplanetary disk; described in detail a new numerical model of this 
equation set; and shown results of vortex formation in small-domain, short-time simulations. 
The baroclinic term is responsible for the vortex production, as it is the only source term 
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in the vorticity equation and the initial vorticity perturbation was zero. Also, increases in 
vorticity and kinetic energy stop when the baroclinic term is turned off. Each vortex has a 
characteristic temperature perturbation around it, due to the temperature advection about 
the vortex; this plays an important role in the baroclinic feedback described in Part II. 

Seventy-four simulations were conducted to test the sensitivity of vortex formation to 
several parameters. It was found that vortex strength increases with: higher background 
temperature (c), stronger background temperature gradients (more negative d), larger initial 
temperature perturbations, higher Reynolds number, and higher resolution. At the highest 
resolution and Reynolds number tested, 256 x 256 and Re = 10^, vortices could be formed 
with a zero initial vorticity perturbation, and an initial temperature perturbation which was 
5% of the background, background temperatures of 200-300K, and background temperature 
gradients that vary radially as r""'^^ and r""'^. 

Additional simulations showed that in this early stage of vortex formation, increasing 
the thermal diffusion weakens vortex formation because the temperature perturbations die 
away too quickly. In Part II we will see that faster thermal diffusion enhances the baroclinic 
feedback and makes vortices grow more quickly, once they have formed. 

The equation set used in this study is anelastic and two-dimensional; this was chosen 
because it is the simplest model which can still include baroclinic effects and radiative cooling 
in a differentially rotating disk. This simple equation set allowed us to create a fast and 
efficient numerical model that was then used to explore a large swath of parameter space 
with hundreds of simulations. By choosing an anelastic 2D model, we limit ourselves to 
dynamics that can occur in that regime, and are blind to processes which may amplify or 
destabilize the vortices that we observe. For example, the anelastic ass umption means that 



the di sk's perturbation velocities are subsonic. The scale analysis in iBarranco fc Marcus 



(120051 ) showed that the horizontal extent of subsonic, compact vortices in a Keplerian shear 
cannot be much greater than the scale height of the disk. Then they performed numerical 
simulations of these compact vortices; when the vortices extend vertically through the disk 
they are unstable to small perturbations and are ripped apart by Keplerian shear. In fact, 
only three-dimensional vortices that are centered above or below the midplane are robust in 
the long-term. 

The vortices in our anelastic simulations are subsonic, compact (typically 0.005-0. lAU), 
and would be a vertical column of vorticity if they were transplanted from their 2D universe 



into a 3D universe. These are like the vortices that iBarranco fc Marcus! (120051 ) found to be 
unstable in 3D. However, their model had a constant background temperature; our vortices 
require a sufficiently steep radial temperature gradient a nd do not develop w i th a constant 



background temperature. So it seems that the vortices in IBarranco fc MarcusI (120051 ) are not 
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due to a baroclinic feedback, even thou gh their model includ e s bar oclinic effects. Regardless 
of the formation process, the work of iBarranco fc MarcusI (120051 ) suggests that columnar 
vortices in a 3D protoplanetary disk are unstable; this provides a strong motivation for 
a follow-up of our study using 3D simulations. It is possible that the parameters that 
influence the baroclinic formation and growth of vortices in 2D (like ours) would affect 
short, stable, off-midplane 3D vortices (like those in lBarranco fc MarcusI l2005l ) in the same 
manner. The current 2D study is important because it elucidates the effects of various 
parameters on the baroclinic instability in an idealized 2D simulation. It also allows one to 
estimate the numerical resolution and effective Reynolds number that would be required in 
a 3D simulation to see the instability. 

By describing the equation set, numerical model, and sensitivity study of parameters, 
this paper lays the groundwork for Part II, which studies the role of the baroclinic feedback in 
growth and longevity of vortices. The simulations in Part II use a quarter annulus domain, 
which is ten times larger than the small domain used here, making it effectively a lower 
resolution model. Thus the findings of the sensitivity analysis apply — in order to form 
vortices, a higher background temperature, steeper background temperature gradient, and 
larger initial perturbation were required. Because of the sensitivity analysis of Part I, it is 
clear that all of these requirements ease with higher resolution and Reynolds number. 

Finally, it is useful to cor npare model pa r amete rs to those inferred from observations 
of real protoplanetary disks. iBeckwith et al.l fll990l ) surveyed 86 pre-main-sequence stars 
in the Taurus-Aurign dark clouds, and found that disk temperatures range from 80-300K 
and vary radially as r'^ with d < —0.5. Our simulations produced vortices with background 
temperatures of 200-300K when d = —0.5, indicating that vortices in our model can be 
generated under realistic conditions. 
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A. Details of the Numerical Model 
A.l. Time-stepping scheme 



The third order, three stage Runge-Kutta scheme by ISpalart et al.l (jl99ll ) is used to 
solve the model equations (iQl fTT]) . The two Jacobian terms are nonlinear and will be treated 
explicitly, while the linear terms on the right will be treated implicitly. 

Define the Jacobians M(^,C) = (9(^,C/So) and A^(*,^) = (9(*,^) /So so that each 
stage of the Runge-Kutta routine is 

(1 + PAtT-^ - /JAtKeV^) = (1 - aAtr-^ + aAtKeV^) 6""-^ 

-7AtA^"-i - 6AtN''-^ (Al) 
(1 - /3Atz/eV') C'" = (1 + aAtz/eV') C'""' 

-7AtAf"-i - 6AtM^-^ + ^— ^ p—- + a—— A2 

r or \ 0(p 0(p J 

where 7,5 are constants of the scheme, n is the current stage, and At is the time step 
(other variables are defined in Sect. [2])- At the beginning of this stage, all n — 1 and n — 2 
variables are known. (At stage 1 the n — 2 variables are zero). Then 6'"^ and C'" are calculated 
by inverting operators of the form (a + feV^) in spectral space (Sect. IA.2I) . The perturbation 
stream function is found by inverting the Laplacian-like operator in iQ in spectral space 
(Sect. IA.3p . Finally, M" and N"' are found by transforming derivatives of 6'^, and 
into physical space, computing the Jacobians, and transforming back into spectral space. 



A. 2. Inverting (a + 6V^) 

Consider the general problem 

(a + 6V')t/(r,0) = /2(r,0), (A3) 

which in polar coordinates is 

a + b— + -— + -— ]u = R (A4) 

The radial basis functions are Chebyshev polynomials, which have a natural domain of 
z G [—1, 1]. To transform from an annular domain with r G [ri,r2], define / = (^2 — Ti)/2 
and g = {ri + r2)/2 so that r = f z + g. By the chain rule drU = dzU //so that (lA4p is now 

Pdz^^ f{fz + g)dz^ {fz + gyd<P')^ " ^^^^ 
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In order to invert this operator in spectral space, expand the variables as 

p K 



(A6) 



p=0 k=0 

where Tp is the p*^ Chebyshev mode. Azimuthal derivatives are simply d"^ = —k"^, so that 
each Fourier mode is independent. The radial derivatives are much more complicated; the 
products of z and d erivatives of z which app ear in flA5D can be represe nted by summations 
with higher modes JCottlieb &: Orszaej [l977l p. 159, Ij.C. MasonI l2002l . p. 32) so that flA5l) 
can be written as Au^ = Brfc where A is an upper triangular matrix, B is a pentadiagonal 
matrix, and and hold the Chebyshev coefficients of U and R for each Fourier mode k. 
Because of the particular form of A, we may add and subtract the rows of Au^ = Br^ in a 
clever way to reduce A to a ba nded diagonal o f width six. Full details and lengthy algebraic 
manipulations can be found in iPeterseru (120041 ). 



A. 3. The stream function 

After the vorticity is found using (lA2p . the stream function is found by inverting 
the Laplacian-like operator in ([9]) in spectral space. This inversion is difficult if the surface 
density So(r) is an arbitrary function of r, but simplified if it is a power function Sq = ar^. 
Then iQ can be rewritten as 

/ 5^ 1 - d d 1 ,„ 

\dr'^ ~^ r dr ~^ r'^dcj)^) ^ ^'^ 

which is nearly the same form as (]A4p . The product So(r)(^'" must be computed in physical 
space and then transformed to spectral space for the inversion. 



A. 4. Boundary conditions - the tau method 



Within spectral methods, there are three m ethods of satisf ying the boundary conditions: 
the pseudospectral, Galerkin, and tau methods (lFornbergj|l996l . p 162). In the pseudospectral 
method, basis functions are chosen that each naturally conform to the boundary conditions. 
In our problem, the azimuthal direction employs Fourier basis functions, which are naturally 
periodic. In the Galerkin method, a new set of basis functions which satisfy the boundary 
condition are created from linear combinations of the old basis functions. 



In the tau method, which we use in the radial, the two highest mode equations in 
Aufc = Brfc are replaced by equations which enforce the boundary conditions. For Dirichlet 
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boundary conditions U{ri,(f)) = C{(f)) and U{r2,(f)) = D{(j)) the corresponding tau lines are 
p p 

^ Un,pTp{-l) = Cfc, ^ Un,pTp{l) = 4, k = O...K (A8) 

p=0 p=0 

and for Neumann boundary conditions drU{ri,(f)) = C{(j)) and drU{r2,(f)) = D{(f)) they are 
p p 

Mn,pi9^Tp(-l) = Cfc, ^ M,i,p(9^Tp(l) = dk, k = 0...K, (A9) 

p=0 p=0 

where Ck and (ifc are the Fourier coefficients of C(0) and D{(f)). In practice C and D are 
constants, so all modes are zero except cq and do- Dirichlet boundary conditions are used for 
impermeable, stress-free boundaries, and Neumann are used for zero thermal flux boundaries. 



A. 5. Stress-free boundary 



The ii npermeable, stress- free radial boundary is enforced using the influence matrix 
technique (jPeyretl |2002| ) . Stress- free boundaries in polar coordinates are complicated by 
the form of the stress tensor so that the conditions for vorticity and stream function are 
coupled. The momentum stream function is defined as EqU = —V x so that the radial 
and azimuthal background velocities are u' = — c?<^\l/'/rEo and v' = —dr'^'/T.Q, and the 
vorticity is defined in fllOl) . An impenetrable boundary implies that u' = 0, which can easily 
be satisfied by requiring that be a constant along the boundary. Stress-free boundaries 
mean that each component of the stress tensor 

1 du' dv' 
r dd) dr 



a 



fx 



2^— 

dr 



V 

r 



Idu' 



dv' 
dr 



V 

r 



2 dv' 2u' 



(AlO) 



is zero (lArisI (119621 ). p. 181). A stream function which is constant along the boundary in 
imphes that c?0\l/ = 0, so that u' = 0, dru' = 0, d^u' = and d^pv' = on the boundary. The 
remaining condition in the stress tensor is 



dv' 
dr 



V 

r 



0. 



(All) 



The stress-free impermeable boundary conditions in terms of stream function and vor- 
ticity are 

^' = 0, (A12) 
2 d^' 



rSn dr 







(A13) 
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where the second condition comes from flAllI) . These coupled equations are more comphcated 
than Dirichlet or Neumann boundary conditions described in the previous section. They 
require the influence matrix technique, where each unknown is written as the sum of a 
particular and homogeneous solutions 

where a and (3 are unknowns to be solved for, superscript p denotes the particular solution 
and superscripts + and — denote the homogeneous solutions. A linear combination of two 
homogeneous solutions is required to construct the desired values at the inner and outer 
boundary of the annulus. To that end, let r~ and be the radial location of the inner 
and outer boundary, and require C~('"~) = 1 C'''('"''') = 1- These are the only functions 
which contribute the the value of ( on the boundary, so set C~('"''') = 0, C^{r~) = 0, and 
(^{r"^) = 0. The stream function should always be zero at the boundary, so require that 
^p(r±) = and *±(r±) = 0. 

Boundary conditions are enforced using the Tau method during the inversion of Laplacian- 
type operators. The vorticity is solved for in each stage of Runge-Kutta time-stepping using 
[2]), which can be written as (a + 6V^)C' = R where R is the right-hand side. The stream 



function is computed from the vorticity using Q, which will be abbreviated here as Vl^' = C- 
The particular solutions must satisfy 

ia + b\/^)C = R, V|*P = C*' (A15) 
while the homogeneous solutions must satisfy 

(a + 6V')C^ = 0, V|^± = 0, (A16) 
so that the full solutions (1A14P satisfy 

{a + bV')C = R, V|vl/' = C'. (A17) 

The only condition which remains to be enforced is the stress-free boundary condition 
(lAlSp . which is what caused all of this trouble in the first place. It can be rewritten in terms 
of the particular and homogeneous solutions flA14p as 



rSo dr J \ rSo dr J rSo dr 

This must be enforced at the inner and outer boundaries, r~ and r"*", resulting in the 2x2 
system of equations with unknowns a and jS. The final vorticity and stream function are 
then constructed with a and (3 as in flA14p . 
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In practice this influence matrix technique is performed for each Fourier mode fc, so that 
f lA14p is now 



*1 



k 



The full procedure conducted at every Runge-Kutta stage is 

1. Solve (a + bV^)C = R and V|^p = for and ^p. 

2. Solve (a + hV'^X^ = and V|^^ = for and 

3. For each Fourier mode k, solve 



2 d^l 



2 

rSo dr 
2 9*+ 



for and Pk- 
4. Reconstruct the full solution using (lA19p . 



Pk 



2 9^ 
2 9^ 



(A19) 



(A20) 
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Fig. 1. — Background temperature Tq (left) and potential temperature (right) profiles 
for simulations TG (top) and T (bottom), where To(r) = c{r /9.5AUY, and d is shown. 
Simulations A, B, and Tau use d = —0.75 and P uses d = —0.5. 
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Fig. 2. — The perturbation vorticity (top) and perturbation temperature T' (bottom) 
from simulation A, with a radial extent of 0.5 AU on a l/64th annulus. The time t refers 
to the orbital period in the middle of the annulus. All simulations begin with zero vortic- 
ity perturbation and a random initial temperature perturbation. Vorticity is immediately 
produced by the baroclinic term, while the initial temperature perturbations shear out and 
decay. 




Fig. 3. — Labels as in Fig. [2], but for later time. After the initial shearing period, thin strips 
of vorticity roll up to form anticyclonic vortices, while the temperature field continues to 
decay. 
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pert Vort, t=2.2 orb per 




Fig. 4. — Close-up view (0.4AU by 0.4 AU) of perturbation variables in simulation A. Four 
anticyclones are labeled (top), with their corresponding temperature perturbations (bottom). 
Each vortex advects warmer fluid inward (towards the left) and colder fluid outward, in a 
clockwise direction. 
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max pert vorticity, rad/yr Kinetic Energy, J/m^ max pert Temperature, K 




0.1 0.5 1 5 0.1 0.5 1 5 0.1 0.5 1 5 



time, orb per time, orb per time, orb per 

Fig. 5. — Comparison of maximum perturbation vorticity \('\ (left), perturbation kinetic 
energy (center), and maximum perturbation temperature \T'\ (right) for simulation series B, 
where the baroclinic term was turned off at the time indicated. When the baroclinic term is 
turned off, there is no longer a source of vorticity. This shows that the vortices are produced 
by the baroclinic instability. The reference simulation (tojj=never) is simulation A. 
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Fig. 6. — Plots showing that vortex strength increases with either initial temperature per- 
turbation, Reynolds number, or resolution. The contour plots show the minimum vorticity 
after five orbital periods, as measured from simulation P with parameter values at each 
asterick. These contours correspond to vortex strength, as shown in the four snapshots of 
perturbation vorticity (' a.t t = 5, where points A to D are labelled on the contour plot. 
As simulations become more realistic with higher resolution and higher Reynolds number, a 
smaller initial perturbation is requried to initiate vortices. 
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Fig. 7. — Plots showing that vortex strength increases with either the radial temperature 
gradient, Reynolds number, or resolution (simulation TG). Plot description is the same as 
Fig. [6l This shows that higher resolution and higher Reynolds numbers allow less negative 
radial temperature gradients to initiate vortices. 
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Fig. 8. — Plots showing that vortex strength increases with the background temperature, 
Reynolds number, or resolution (simulation T). Plot description is the same as Fig. [61 
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Fig. 9. — Data for simulation set Tau, where the radiative coohng rate is varied between 
r = 1 and r = 100 orbital periods. With r = 1 the disk cools faster (right), which weakens 
temperature gradients so that less vorticity is produced by the baroclinic feedback. During 
vortex formation, vortices are strongest as r — oo, i.e., in the limit of no radiative cooling. 
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Fig. 10. — Data for simulation set Pe. When thermal diffusion is higher (lower Peclet number, 
Pe = le4) the disk cools faster (right). This weakens temperature gradients, making the 
baroclinic feedback weaker, so that growth in vorticity (left) and perturbation kinetic energy 
(center) is much slower. 
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name 


d 


c 


init. pert. 


r 


grid 


Re 


Pe 


A 


-0.75 


300 


0.05 


100 


256^ 


108 


10^ 


B 


-0.75 


300 


0.05 


100 


256^ 


108 


10^ 


TG 


-0.25- -2 


300 


0.05 


1 


128^, 256^ 


10^-10^ 


10^ 


T 


-0.5 


200-400 


0.05 


1 


128^, 256^ 


10^-10^ 


10^ 


P 


-0.5 


300 


0.02-0.1 


1 


128^, 256^ 


10^-10^ 


10^ 


Tau 


-0.75 


300 


0.05 


1-100 


256^ 


108 


10^ 


Pe 


-0.75 


300 


0.05 


100 


256^ 


108 


lo^lo^ 



Table 1: Model parameters for the numerical simulations discussed in this paper. Here d is 
the power on the background temperature function, c is the background temperature at the 
inner radius, init. pert, is the magnitude of the initial temperature perturbation relative to 
the background, r is the radiative cooling time. Re is the Reynolds number, and Pe is the 
Peclet number. Simulation B is identical to A except that the baroclinic term is turned off 
during the simulation. 



